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We study the process of vortex nucleation in rotating two-dimensional BEC confined in a harmonic 
trap. We show that, within the Gross-Pitaevskii theory with the boundary condition of vanishing 
of the order parameter at infinity, topological defects nucleation occurs via the creation of vortex- 
antivortex pairs far from the cloud center, where the modulus of the order parameter is small. Then, 
vortices move towards the center of the cloud and antivortices move in the opposite direction but 
never disappear. We also discuss the role of surface modes in this process. 
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I. INTRODUCTION 
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Quantized vortices in Bose-Einstein condensates (BEC) of alkali atom gases attract a considerable current interest. 
One of the most interesting problems of vort ex phys ics is the nucleation of top olog ical defects. A lot of experimental 
H H H || and theoretical |E S S S, 0, E3, IHl E3, E3, EH [13, EH Ea E3, HH works deal with this subject. 
^ ■ Phases with different number of vortices are usually separated by large energetic barriers. Therefore, the transition 
to the vortex state occurs at higher rotation frequency than the thermodynamical one determined by the balance of 
energies of vortex- free and one- vortex states. The same is valid for the well-studied case of superconductors [23ll2^| . 
However, BEC clouds have no well-defined surface and, therefore, the process of vortex formation must be different 
I i from that in superconductors, where vortex lines penetrate from the sample surface. Within the Gross-Pitaevskii 
^ ' theory, the boundary condition for the order parameter in the harmonically trapped BEC reads 

V>(r^oo) = 0. (1) 

In the lil of strongly interacting gas, a Thomas-Fermi approximation is often used. In this approach, the Thomas- 
Fermi boundary, where tp vanishes identically, has to be assumed, and vortices are usually treated as point-like 
, objects. Under these assumptions, one can calculate the critical rotation frequency, which suppresses the surface 
I/"") • barrier preven ting vortex penetration through the Thomas- Fermi boundary (jj, llll Il2| . Another approach to the 
problem [H [tj 

M m is based on the analysis of stability of the anally- S yru7n~icTorte X -hee state with respect 
^sO • to the perturbations of the order parameter tp, and it was shown that the nonvortex state becomes unstable with 
respect to surface modes, which leads to the well-known Landau instability. It is believed that this instability results 
in the nucleation of vortices. The third approach to the pro blem of vortex formation is by numerical solutions of 
the time-dependent Gross-Pitaevskii (GP) equation |ol Il4 fl5l IT^. supplemented by the boundary condition of 
varnishing of ip at some boundary. 

The aim of the present paper is to study a vortex nucleation process taking into account consistently boundary 
condition (1) for the order parameter. We consider the case, when the system is spined up quasistatically and 
adiabatically. No boundary, which can serve as a gate for vortex penetration is assumed. We use a semi-analytical 
model [2j|, which is applicable for the lil of diluted gas. We show that nucleation of vortices in the model with 
boundary condition (1) occurs via the creation of vortex-antivortex pairs far from the cloud center, in regions, where 
\tp\ is small. Vortices then move toward the cloud center to decrease the energy of the system, whereas antivortices 
go to the opposite direction. This mechanism for the vortex nucleation is natural for the case of a system, which has 
no boundary and where the total topological charge has to be conserved. Also studied is the role of surface modes 
in the process of vortex nucleation. We show that the critical velocity for vortex formation is slightly higher than 
the Landau critical velocity, and this is in a agreement with Anglin's results [Til IT^| obtained by the very different 
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method. 



II. MODEL 

We consider the case of quasi two-dimensional condensate. The dimensionless GP energy functional is given by: 



F = I rdr 



d^|W| 2 + y M 2 + ^H 4 + i»P jg), (2) 
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where the integration is performed over the area of the system, (r, ip) are polar coordinates, N is the number of atoms, 
lu is the rotation frequency, p = an z is the gas parameter, a and n z are scattering length and concentration of particles 
along the cloud axis, respectively. Distances and rotation frequencies are measured in units of the oscillator length and 
the trapping frequency, respectively. The normalization condition for the order parameter reads J rdr J dip \ip\ = N. 



In general case, ip can be represented as a Fourier expansion 



(3) 



In the lil of noninteracting gas (p — 0), it can be easily shown by substitution of Eq. (3) to the GP equation 
that each function fi(r) coincides with the eigen function for the harmonic oscillator corresponding to the angular 

momentum I. These functions have the Gaussian profile ~ r l exp ^— ^-J . Therefore, one can assume that this Gaussian 

approximation remains accurate in the case of weakly interacting dilute gas. The accuracy can be improved if we 
introduce a variational parameter R[ characterizing the spatial extend of fi(r). Finally, our ansatz for fi(r) has a 
form: 



exp 



2Rf 



~ 1><PI 



(4) 



where C/, i?/, and (pi can be found from the condition of minimum of the energy (2), Ci is a real number. This 
approach was used for the first time in Ref. j2^, see also Ref. [2(||22j]. In the present paper, we restrict ourselves to 
the interval < p < 10, since, according to our estimates, the Gaussian approximation becomes inaccurate at higher 
p. To verify this, we compared known numerical results for the thermodynamical critical rotation frequency with 
those obtained by our method. 

Now we substitute Eqs. (3) and (4) to Eq. (2) and after the integration we obtain: 



where 



l l l>k 

+4 hkk m ClClC m 8i +m ^k COS (</>; + (p m - 2<Pk) 

l>k>m 

+8 ^ ' IlkmnClCkC m C n 5i^k, 

m+n COS ((pi + (f) k - (f) m - (j) n ) , 



l>k>m>n 



at = -T(l + 2) (1 + Rf) + nR?T(l + FM, 



(5) 



(6) 



2ir 2 p^ I + m + n + k 
N ^ 2 
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Rn 



Rlkmn = V% (Ri 2 + Rk 2 + Rm 



R. 



T(l) is a gamma function. Normalization condition is now given by 

nJ2C?R 2 T(l + l) = N. 



Rlkr 



Rn 



(7) 
(8) 
(9) 



Values of parameters i?/, Ci and (pi can be found from the minimum of the energy (5) taking into account Eq. (9). 

For instance, for the axially-symmetric vortex-free state, Co = ^JN/ttRq, Rq = (1 + 2p) 1 ^ 4 , and C\ = at I ^ 1. 
In the vicinity of the local minimum, the deviation of the energy from the minimum can be represented as 

AF = Y / Ai k AC l AC k , (10) 

Ik 

where AC; is a deviation of Ci from its equilibrium value. In Eq. (10) we kept the terms up to the second order of 
AC;. The stability criterion of the solution is a positiveness of all main minors of the matrix 
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III. STABILITY OF NONVORTEX STATE 



First, we study the stability of the vortex- free state and after this we will go beyond this linear approach. We 
consider the situation, when the system is spined up quasistatically. For the axially-symmetric nonvortex state, the 
matrix \\Aik\\ is diagonal. In this case, it is easy to find u>, at which vortex-free phase becomes unstable with respect 
to the surface mode with the angular momentum I: 



IT (I 
1 

4 



1 + RV 



(z - i)r(/ + 1) 

T(l + 2) + 



4(1 + 2p)'/ 2 r<7 + 1) 



(11) 



Value of Ri can be found from the condition of minimum of u)i ns t(l,p) with respect to Ri. In Fig. 1 (a) we plot the 
dependences of u>mst{l,p) on p for different I. The resulting function Wi ns t{p) is determined by the minimal value of 
Uinstihp) with respect to the quantized I. This function is shown in Fig. 1 (b) and represents a well-known Landau 
criterion for the angular- momentum transfer to the superfluid system. From Fig. 1 (b) one can see that value of li ns t, 
which leads to the instability, increases with increase of gas parameter p. Vortex-free phase is always stable with 
respect to the generation of harmonic 1=1. At low p values, the instability with respect to the harmonic I = 3 occurs 
at lower ui than that for 1 — 2. In Fig. 1 (b), we also presented the thermodynamical critical rotation frequency as a 
function of p, which we calculated by the comparison of energies of vortex-free and single-vortex states. It is much 
lower than oJinst, as can be expected. 



IV. NUCLEATION OF VORTICES 



Now we go beyond the linear analysis. It follows from the results of the previous Section that, in the instability 
point, an admixture of a harmonic with I = li ns t appears in the order parameter expansion (3). In fact, the order 
parameter must contain also all other harmonics with Z's divisible by h ns t, which are induced by the component with 
I = linst, as can be easily seen from the GP equation. However, the main contribution to the energy is given by 
just two first terms with I = and I — Unst- Therefore, we will take into account only these two harmonics. The 
component with I = li ns t gives ^ nS f-fold modulation of the modulus of the order parameter in the azimuthal direction 
and it is responsible for the creation of vortices. We calculate the energy and tp in the two-harmonic approximation at 
values of w slightly exceeding w inst and study their evolution with changing of u>. At each u> we also check the stability 
of the obtained solution with respect to the nucleation of other harmonics with various Z's, not divisible by li ns t- The 
two-harmonic approximation was used widely for the analysis of vortex states in mesoscopic superconducting discs 
within the Ginzburg-Landau theory, see e.g. Refs. jH, E^, EH EH ■ Note that we also tried to take into account 
several additional components of tj) with Vs divisible by hnst (2hnst, 3h ns t, ^hnst), but this does not change our results 
significantly thus verifying good accuracy of the two-harmonic approximation in our case. 

First, we discuss the simplest case of Unst = 3. In the instability point, an admixture of harmonic with I = 3 appears. 
Surface mode leads to the three-fold modulation of |?/>| in the azimuthal direction. In Fig. 2 we show the evolution of 
\tp\ far from the system center along one of three directions, where \tp\ is minimal, at p = 1. In the beginning, |^| in 
this direction is a monotonic function of r. With further increasing of u>, a local minimum appears at the dependence 
of ip( r )i see Fig- 2 (a). This minimum becomes deeper and deeper and then vanishes at some point, as shown on 
Fig. 2 (b). If we proceed with spinning up the system, the minimum splits to the vortex-antivortex pair, see Fig. 2 
(c). If we go around the center of the vortex, the phase of the order parameter changes by 2-7T, and around antivortex 
it changes by — 2ir. Thus, we can conclude that these zeros of the order parameter indeed correspond to the vortex 
and antivortex. Increasing of u> leads to the motion of vortex toward the system center and antivortex in the opposite 
direction, where \tp\ is very small. The same is happening along the other two directions. As a result, three vortices 
penetrate the inner part of a cloud in a symmetrical pattern, whereas three antivortices move far from the cloud 
center but never disappear. The spot, where vortex-antivortex pair nucleates, is characterized by a small value of \tp\, 
approximately, fifty times smaller than \ip\ in the cloud center. Therefore, the formation of pairs is invisible on the 
large length-scale. Also, the interval of to between the point of the instability and the vortex state is very narrow, of 
the order of 10 -3 of trapping frequency. This interval is occupied by the vortex-free state with the order parameter 
modulated in the azimuthal direction. We believe that vortex-antivortex pairs are not an artifact of our method, since 
they appear even in the lil p << 1, where Gaussian approximation is very accurate. Also, this mechanism of vortex 
nucleation is rather natural for the continuos system without any borders, as reflected by the boundary condition (1) 
for the order parameter. 
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The similar scenario for the vortex nucleation is realized at values of p corresponding to Unst — 5. In this case, 
five vortices penetrate the inner part of the cloud and five antivortices move in the opposite direction. Thus, surface 
modes induce ripples of the order parameter and then ripples are naturally transformed to the same number of vortex- 
antivortex pairs. However, when Unst is even, li ns t — 4 or 6, the scenario is different. At first, \tp\ becomes modulated 
in the azimuthal direction and with increasing of to, the minima of "0 become more and more pronounced exactly as 
at odd Unst- But just before the formation of Unst vortices this phase suddenly becomes unstable with respect to the 
nucleation of harmonic with I — Unst/^- The stationary state after this instability corresponds to the local minimum 
of energy with tj), which is given by a superposition of harmonics with I — 0, ^„ s t/2, and Unst- At 1 < p < 1.5 
(knst = 4), this phase is also unstable with respect to the nucleation of the harmonic with 1 = 1. This implies that 
the final state contains one vortex in the center of the cloud and one antivortex at infinity. Thus, according to our 
model, for even values of U ns t, not all ripples of ip are transformed to vortex-antivortex pairs. Number of these pairs 
is equal to 3 for Unst — 6 and 1 or 2 for Unst = 4 depending on p. The role of surface modes in this case is to facilitate 
the penetration of one or more number of vortices to the inner part of the cloud. The similar mechanism for vortex 
penetration was analyzed in Ref. Q for the particular case of a quadrupole mode, / = 2, which facilitates penetration 
of a single vortex. No antivortices were found in this case, since the Thomas-Fermi boundary was assumed. Our 
result is more general. It shows that penetration of vortices can be related to the induction of different surface 
modes depending on p. After the cascade of various harmonics generation, multiple- vortex penetration can occur, not 
necessarily one vortex enters the cloud. In our model, the barrier preventing the separation of vortex and antivortex 
(or surface barrier in models assuming Thomas-Fermi boundary) is destroyed at critical frequencies sli ghtl y exceeding 
the Landau critical frequency for the surface modes. Thus, our model supports the Anglin's results [Til Il2|. which 
were obtained by using completely different method, namely, the hydrodynamic approach taking into account spatial 
variation of the order parameter in the Thomas-Fermi surface layer. 

The mechanism for the vortex nucleation via vortex-antivortex pairs is a specific feature of infinite systems. One 
can see here some analogy with the BKT transition, where topological defects also nucleate by pairs, and the total 
topological charge is always zero. In finite systems, topological defects usually penetrate through the surface, but 
here no surface is assumed, since condition (1) should be fulfilled. Therefore, a natural way to create vortices is 
via vortex-antivortex pairs, after which antivortices are pushed to the regions with small \tp\, since their presence in 
regions with high is energetically unfavorable. Therefore, antivortices are practically undetectable experimentally. 
An important consequence of this result is that the topological charge of BEC with boundary condition (1) is always 
equal to zero. This means that for any stable configuration of vortices in the inner part of the cloud there is a 
configuration of antivortices at the periphery of the system, where the d2ity of particles is very low. Indeed, using our 
method we calculated ip(r, (p) for different numbers of vortices and we always found the same number of antivortices 
in the system. Note that some additional topological defects in regions of small were reported in Refs. 0,0>E1j 
where int2ive numerical methods were used for the solution of the GP equation. They were called "ghost vortices" . 
The correspondence between the antivortices in our model and "ghost vortices" is an open question. 

It can be also predicted that in 3D case vortex loops nucleate in rotating BEC instead of vortex-antivortex pairs. 
With increasing ui, lateral dim2ions of loops grow up and finally 'vortex' parts of loops penetrate the inner part of 
cloud and 'antivortex' parts go to the opposite direction. This effect is related to the bending of vortex lines found in 
3D numerical simulations ^3, . It is also interesting to analyze the stability of knotted GP solutions in this case. 
In 3D superconductors, vortices also penetrate through the surface not as straight lines, but as half-loops pil l32|. 
Physically, this is due to the fact that the vortex energy is proportional to its length and therefore a creation of vortex 
from the point nucleus is easier than from the line one. BKT transition in 3D case also occurs via the generation of 
vortex loops |3^| , see also [34j . 

Now we briefly discuss the case of finite cloud. We perform the same calculations, but using functions 



fi(r) = a 



+ A 



1+2 



exp 



2Rf 



(12) 



where (3i is an additional parameter, which is chosen in such a way as to meet the boundary condition fi(Rb) = 0, is 
the radius of the system. We found that topological defects nucleate via vortex-antivortex pairs, if Rb is approximately 
two times larger than the Thomas-Fermi radius. With further spinning up of the system, some of antivortices leave 
the system and the topological charge is not necessarily zero. When R^ is less, vortices penetrate via the surface as 
usual, and no antivortices appear. 
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V. CONCLUSIONS 



We studied a nucleation of topological defects in rotating BEC in 2D case within the Gross-Pitacvskii theory with 
boundary condition of varnishing of the order parameter in infinity. We found that in this model, with increasing 
slowly of the rotation frequency, vortex-antivortex pairs appear in the system. Vortices move to the inner part of 
the cloud and antivortices are pushed to the opposite direction. A number of pairs depends on the gas parameter. 
Topological charge of the system is always equal to zero, and for any configuration of vortices there is a configuration of 
antivortices far from the cloud center. For 3D case, we predict that vortex loops nucleate instead of vortex-antivortex 
pairs. Also studied is the role of surface modes in the process of vortex formation. We revealed two scenarios: either 
all the ripples of the order parameter induced by surface modes are transformed to the vortices or surface modes 
facilitate penetration of vortices through the cascade of various harmonics generation. A multiple penetration of 
vortices is possible. In the finite-size clouds, vortices can either penetrate through the surface or nucleate via vortex- 
antivortex pairs depending on the system size. Vortex nucleation occurs at critical velocities slightly exceeding the 
Landau critical velocity for the surface modes. 
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VI. FIGURES 

Fig. 1. The dependences of rotation frequences LUi ns t(l), at which the axially-symmetric nonvortex state becomes 
unstable, on gas parameter p at different values of quantized angular momentum I (Fig. 1 (a)). Fig. 1 (b) shows the 
resulting function, which is given by minimum of uj inst (l) with respect to /. Dot lines on Fig. 1 (b) are giudes for eyes 
indicating the boundaries between regions corresponding to different values of k ns t ■ Dashed line on Fig. 1 (b) is the 
thermodynamical critical frequency. 

Fig. 2. Evolution of the modulus of the order parameter in the vicinity of the spot, where vortex-antivortex pair 
nucleates. Fig. 2 (a) corresponds to the rotation speed 0.77405, Fig. 2 (b) to 0.77417, and Fig. 2 (c) to 0.77419. 
Number of particles is 10000, p = 1. 
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